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Abstract — A new method of matrix spectral factorization is proposed 
which reliably computes an approximate spectral factor of any matrix spec- 
tral density that admits spectral factorization. 

Index Terms — Matrix spectral factorization algorithm. 



1. Introduction 

Spectral factorization plays a prominent role in a wide range of fields in Communi- 
cations, System Theory, Control Engineering and so on. In the scalar case arising for 
single input and single output systems, the factorization problem is relatively easy and 
several classical algorithms exist to tackle it (see the survey paper [17]) together with 
reliable information on their software implementations [8]. There are also some re- 
cent claims as to their improvement [2] . Matrix spectral factorization which arises for 
multi-dimensional systems is essentially more difficult (see Sect. 2, where the mathe- 
matical reasons of this fact are explained). Since Wiener's original efforts [19] to create 
a sound computational method of such factorization, tens of different algorithms have 
appeared in the literature (see the survey papers [16], [17] and the references therein), 
but none of them is thought to have an essential superiority over all others (see [16, 
p. 1077], [14, p. 206]). Besides, most of these algorithms impose extra restrictions 
on matrix spectral densities (e.g., to be real or rational or nonsingular on the bound- 
ary), while the Paley- Wiener necessary and sufficient condition (see (2)) will do for 
the existence of spectral factorization (see Sect. 2). 

In the present paper, a new computational method of matrix spectral factorization 
is developed. The proposed algorithm can be applied to any matrix spectral density 
satisfying the Paley- Wiener condition. It should be said that the branch of mathemat- 
ics where the spectral factorization problem is posed in its general non-rational setting 
(see Sect. 2) is the theory of Hardy spaces (see Sect. 3), and this method is completely 
worked out in the framework of Hardy spaces, which added to its effectiveness. 

To describe our method of r x r matrix spectral factorization in a few words, it 
carries out spectral factorization of m x m left-upper submatrices step-by-step, m = 
l,2,...,r. It is shown that in this process the decisive role is played by unitary 
matrix functions of certain structure (see Theorem 1), which removes many technical 
difficulties connected with computation. The explicit construction of such matrices 
in Theorem 2 is an essential component of the algorithm. A close relationship of 
these unitary matrix functions with compactly supported wavelets has recently been 



^This paper includes the detailed proofs for an innovative method for matrix spectral factorization 
that can be used in numerous applications, including Filtering, Data Compression, and Wireless 
Communications. A U.S. patent application has been submitted for this innovation through the 
Technology Commercialization Center of the University of Maryland. 
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discovered, which makes it possible to construct compact wavelets in a fast and reliable 
way and to completely parameterize them (see [6]). 

Preliminary numerical simulations confirm the potential of the proposed algorithm 
(see Sect. 7). 

The algorithm was announced in [3] and, for second order matrices, described in 
[12]. 

2. Formulation of the problem 

A series of papers [18], [19], [10], [11] led to the following 
Wiener Matrix Spectral Factorization Theorem: Let 



(1) 


S{t) = 


S2l(t) 


Sl2{t) ■■ 
S22it) ■■ 


■ Sir{t)\ 

■ S2r{t) 






\Srlit) 


Sr2{t) ■ ■ 





t G T, be a positive definite (a.e.) integrable matrix function, < S{t) e Li(T), which 
satisfies the condition 

(2) logdet5(t) e Li(T). 
Then it admits a spectral factorization 

(3) S{t)^S^{t)S-{t) = S-'it){S^t)y, 

where is an r x r outer analytic matrix function from the Hardy space H2 and 
S~{z) = (^S~^{l/z)^ , \z\ > 1. It is assumed that (3) holds a.e. on T. (The factorization 

(3) is called left since the analytic inside T factor stands on the left-hand side. The 
right spectral factorization of S can be obtained by the left factorization of 5"^.) 

The sufficient condition (2) is also a necessary one for the factorization (3) to exist 
(see Sect. 3). 

A spectral factor S~^{z) is unique up to a constant right unitary multiplier (see, 
e.g., [5]), and the unique spectral factor with an additional requirement that 5'+(0) be 
positive definite is called canonical. 

After the proof of the existence of matrix spectral factorization, the computation 
of the spectral factor for a given matrix spectral density has become a challenging 
problem due to its applications in practice. 

In the scalar case, r — 1, the canonical spectral factor 5"+ e H2 can be explicitly 
written by the formula (see, e.g., [20; VII, 7.33]) 

/ 1 /•27r „id , \ 

(4) 5+(.) = exp -^JogS{e^')d9j 

and it is relatively easy to compute approximately. However, there is no ana- 
log of this formula in the matrix case because, generally speaking, e"^"*"^ 7^ e^e^ for 
non-commutative matrices A and B. This is the main reason for which the approx- 
imate computation of the spectral factor S'^ in (3) for the matrix spectral density 
(1) is essentially more difficult. The present paper provides an algorithm for such 
computation. 
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The proposed method does not contribute to the improvement of (numerical) scalar 
spectral factorization, but employs it to fulfill matrix spectral factorization. 

3. Notation and Conventions 

Let D = {2; e C : 1^1 < 1}, and T = c?D be the unit circle. As usual, Lp = I/p(T), 
< p < 00, denotes the Lebesgue space of p-integrable complex functions defined on 
T. Hp — Hp(B), < p < 00, is the Hardy space of analytic functions in D , 

Hp^lfe A(D) : sup / \f{re''^)\Pde < 00] 

I r<l Jo ) 

{Hoo is the space of bounded analytic functions), and = L^C^) denotes the class 
of their boundary functions. (All the relations for functions from L^{T) or Lp(T) 
are assumed to hold almost everywhere.) Since there is a one-to-one correspondence 
between Hp{3) and (T), p > (see, e.g., [20; Vll, 7.25]), we naturally regard these 
two classes as identical, and thus we can speak about the values of / G L^{T) inside 
the unit circle. Furthermore, we always use the argument t for the functions defined 
on T and the argument z for the functions defined in D, so that the boundary function 
of / = f{z) e Hp{B) is denoted by / = f{t) G L+(T) and we write f{z)U=t = f{t) 
when we wish to point out this fact. If we write only /, its domain will be clear from 
the context. 

We have \og\f{t)\ G Li(T) for each ^ / G i/p, p > (see, e.g., [20; VII, 7.25]), 
which readily implies the necessity of the condition (2) for the factorization (3) to exist 
since Li(T) 3 log | det 5+(t)| = | logdct 5(t). 

The nth Fourier coefficient of an integrable function / G Li(T) is denoted by c„(/). 
For p > 1, I/p (T) coincides with the class of functions from Lp{T) whose Fourier 
coefficients with negative indices are equal to zero. We also deal with -L~(T) = {/ : 
/ G L+(T)} = {/ G Lp{T) : c„(/) = whenever n > 0}. The set of trigonometric 
polynomials is denoted by V, i.e. / G P if / has only a finite number of nonzero Fourier 
coefficients. Also let := V n L^, Vn ■= {f ^ V : Cn{f) = whenever \n\ > N}, 
and P± = n Obviously, f eV^^J eV^. 

For fit) = Er=-oo Cn^" e H^), let P+f{t), P-f{t), and P^/(i) be the projections 
E^o^n^", E°=-ooCn^", and YZ=-NCnt'^^ respectively, on L+(T), L^{T), and Vn- 

The superscript "+" (resp. "— ") of a function /"^ (resp. /~) emphasizes that this 
function belongs to Lp (resp. L'). 

The norms || ■ and || ■ \\h^ are defined in a usual way. 

If M is a matrix, then M denotes the matrix with conjugate entries and M* := . 
If M is positive definite, M > 0, then the unique Mq > that satisfies MqMq — M is 
denoted by -s/M. 

If M is an r X r matrix and m < r, then {M)mxm is assumed be the mxm upper-left 
submatrix of M. 

An r X r matrix U is called unitary if UU* — U*U — I^., where stands for the 
r-dimensional unit matrix. Obviously the entries of a unitary matrix are bounded by 
1. 
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A matrix function M{t) defined on T is called positive definite or unitary if it is 
such for almost all t G T. M{t) is said to belong to some class, say, L^{T) (we write 
M{t) G L+(T)) if its entries belong to this class. Pj^M^t) denotes the matrix function 
whose entries are the projections of the entries of M{t) on Vn- A sequence of matrix 
functions is said to be convergent in Lp-norm if their entries are convergent in this 
norm. 

The class of r x r unitary matrix functions U{t), 

(5) U{t)U*{t) = Ir a.e., 

is denoted by W,.(T), and SUr{T) stands for the subclass of those U{t) G Ur{T) the 
determinants of which are equal to 1, 

(6) det U{t) = 1 a.e. 

The set of outer analytic functions from the Hardy space Hp, p > 0, is denoted by 
Op. RecaU that / G Op if and only if ^ / G -f/p and 

/(-) = ^ • log \fi^'') I d^) ' 1^1 = 1- 

(From this definition and Holder's inequality it follows that if / G Op and g & Og, 
then fg G 0(p+g)/pg .) Clearly, f{z) ^ for each z G D and \f(t)\ > for a. a. t G T if 
/ G Op. The set of functions f & Op which are positive at the origin (which happens 
when c = 1 in the above definition) is denoted by Op. We say that a r x r matrix 
function M(t) G Hp, p > 0, is outer if its determinant belongs to Op/r. This definition 
coincides with some other equivalent definitions of outer matrix functions (see, e.g., 
[11]). M{t) G means that M(0) > in addition. 
f^^f means that converges to / in measure. 

{■,-)m and II ■ llc" denote the usual scalar product and the norm, respectively, in the 
m-dimensional complex space C"*. 
5ij stands for the Kronecker delta. 

To conclude the section, we formulate a simple statement from the Lebesgue integral 
theory in the best suitable form for further references. 

Statement 1. Let /„(t) g UiT), n = 1, 2, . . ||/„(t) - /(t)|U, ^ 0, M„(t) G L^{T), 
Un{t) < I, n = 1,2,..., and Un{t) =^ u{t). Then ||/„(t)M„(t) - /(^)m(^)||l2 ^ 
{see,e.g., [9;§26,T/i. 3]). 



4. Mathematical Background of the Method 

In this section we formulate some statements needed to describe our method. Most 
of the proofs are given in the next sections. 

The uniqueness of spectral factorization (3) mentioned in Sect. 2 means that S~^{z) ■ 
U is also a spectral factor for any (constant) unitary matrix U, and if Si{z) and 5*^(2;) 
are two spectral factors, then Si{z) = S2{z)U for some unitary matrix U (see, e.g., 
[5]). Since for any r x r non-singular matrix S there exists a unique unitary matrix 
U which makes the product SU positive definite (see, e.g. [7; IX §14]), the canonical 
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spectral factor S^{z) (with an additional requirement that 'S'+(0) be positive definite) 
is unique. Namely, 

(7) S^z) = S^{z){S+{0))-',/S+{0){S+{0)y 

for any spectral factor S~^{z). (Other uniqueness restrictions on can be imposed 
so that 5*^(0) would be, for example, lower triangular with positive entries on the 
diagonal.) The following lemma can be applied for the approximation of the canonical 
spectral factor after the approximate computation of an arbitrary spectral factor. 

Lemma 1. Let S~^{t) be a spectral factor of (1) and let 5'+(t) e H2, n = 1,2, . . ., be 
such that 

(8) \\S+{t) - S+{t)\\H, ^0 as 00. 
Then 

(9) \\S+iz){S+iO))-'^S+mS+{0)y - St{z)\\^,^^0 as 00. 

Proof. Since 5*^(0) is non-singular and (8) implies that 5*^(0) 5*^(0), we have 
^5+(0)(^+(0))* - ^5+(0)(!s+(0))* and (^+(0))-' ^ (^+(0))'^ Therefore (9) 
follows from (8) and (7). □ 

The following lemma is used several times throughout the paper. 
Lemma 2. Let M{t) be any m x m matrix function from I/2(T) satisfying 

(10) det M{t) e 02/m C Hy^ = L+JT). 
IfU{t) e SUmi^) is such that 

(11) M{t)U{t) e L+(T) 
holds, then M{t)U{t) is a spectral factor of M{t)M*{t). 
Proof. Taking into account (5), we have 

M{t)U{t) ■ {M{t)U{t)y = M{t)U{t)U*{t)M*{t) = M{t)M*{t). 

In view of (11), MU can be extended inside T. Hence, by virtue of (6), 

dei{MU){z)\^=t = det {M{t)U{t)) = detM(i) = detM{z%^t 

Thus det(MC/)(2;) is an outer analytic function (see (10)) and lemma holds. □ 

This proof gives rise to 

Corollary 1. LetM{t) beanymxm matrix function from L2{T) satisfying det M{t) e 
^f^i^) ^ «5^m(T) is such that (11) holds, then 

(12) det {MU){z) = {det M){z), \z\ < 1. 

It should be pointed out that on the left-hand side of (12) we first extend MU inside 
T and then compute its determinant, while on the right-hand side we first take the 
determinant of M{t) and then extend it inside T. 

The following two theorems play a decisive role in our method. 



Theorem 1. For every m x m matrix function F{t) e L2{T) of the form 



(13) 



( 1 








1 








1 









\ 






1 





vci(^) m m 

where 

(14) 0(t) e L2(T), J = 1, 2, . . . , m - 1, and f+{t) e Ol C L+(T), 

there exists a unique Upit) e 5Wm.(T) o/ the form 



(15) 
(16) 

such that 
(17) 



""m-ial^) '^m-l,2(^) 



'^rn-l,m-li^) '^m-l,mi^) 



Tn,m- 



-lit) 



U 



mm 



it) J 



umeL^iT), i,j ^1,2,. ..,m, 



Fit)U^{t)^F:{t)eO',GLtin 



where F^(t) is the canonical spectral factor of F(t)F*(t) . 

The proof of Theorem 1 relying on the existence of spectral factorization is relatively 
easy (see [3]). The core of the proposed matrix spectral factorization method is the 
constructive proof of Theorem 1 based on the following idea. Wc approximate F{t) in 
L2 cutting off the tails of Fourier expansions of the functions Q{t), j = 1, 2, . . . , m — 1, 
and f~^{t). Namely, for a matrix function of the form (13), (14), let F^^\t) be PArF(t), 
i.e. 



(1^ 



where 



(19) 



/ 1 








\ 


1 














1 



















(At) ■■■ 






N 




N 




= ^c„(/+)r. 


and Cf\t) 






n=0 




n=-N 




is obvious that, 
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or, equivalently, 

(20) ||FW(i) - F{t)\\L, ^0 as N^oo. 

We will multiply (18) by the polynomial unitary matrix function which eliminates 
the Fourier coefficients with negative indices of the product. Furthermore, we prove 
the following theorem for matrix-functions F^^\t), N — 1,2, . . ., which involves the 
limiting case too. 

Theorem 2. (a) Let N be any positive integer, and let a matrix function F^^\t) e Vn 
of the form (18) be such that 

(21) Cf\t) eVN, j^l,2,...,m-l, and /+ ^(i) e V^, /+ )(0) > 0. 
Then there exists and one can explicitly construct 

(22) UpiN){t) e SUm{T) 
of the form (15) such that 

(23) umeV^, i,j = l,2,...,m. 



(24) FW(i)C/pw(i)eP+, 
and 

(25) F^^^Up(N){0) > 0. 

(b) Given an arbitrary sequence of matrix functions F^^\t), N = 1,2,. . ., of the 
form (18), (21) which converges in L2 to F{t) {i.e. (20) holds) of the form (13), (14), 
we have 

(26) \\F^''\t)Up,N,{t) - F+{t)\\H, ^0 as 00. 

Furthermore, the sequence Up(N){t), N = 1,2,..., is convergent in measure. The 
limiting matrix function U{t) G SUmi^) satisfies the conditions imposed on Up{t) in 
Theorem 1, and therefore U{t) coincides with Up{t). Consequently, we have 

(27) f/^,^)(t) ^f/^(t). 

The constructive proof of Theorem 2 (a) given in Sect 5, which computes explicitly 
and in a fast reliable way the coefficients of the functions M^(t) G in (23), is 
the essence of the proposed algorithm. The part (b) of the theorem involves the 
algorithm convergence properties and is proved in Appendix A. We point out the fact 
that Theorem 2 (b) includes also the proof of Theorem 1. 



5. A Constructive Proof of Theorem 2 (a) 

Throughout this section it is assumed that N is fixed and Q(t) :— (j^\t), f'^(t) :— 
/+)(t),and F{t) :=FW(t). 



For given functions Cjif)i 3 = 1, 2, . . . , m — 1, and f~^(t) satisfying (21), we consider 
the following system of m conditions, which plays a key role in the proof, 



(28) 



Ci(t)a;+(t)-/+(t)a:+(t)eP+, 
C2(t)x+(t)-/+(t)x+(t)eP+, 



C™-i(t)x+(t) - /+(t)x+_i(t) e V+, 

^Ci{t)xt{t) + C2(t)4W + • • • + Cm-i{t)xi_,{t) + f+{t)x+{t) e r+, 



where the vector function (^x'l (t) , X2 {t) , 
We say that a vector function 



(i)) is unknown. 
{ut{t),ut{t),...,ui{t)f 



(29) U (t) — y^u,^ i^t;, u,2 v;, • • • , "mW; ^ ''AT 

is a solution of (28) if and only if all the conditions in (28) are satisfied whenever 
xf{t) = u^{t), i = 1,2, ... ,m. Observe that the set of solutions of (28) is a linear 
subspace of m-dimensional vector- valued functions defined on T. 

For the vector function (29), we define the modified vector function u+{t) as 

(30) i^it) = {utit),uUt),...,^)f. 

It is assumed that the modification of (30) is (29). 
We make essential use of the following 

Lemma 3. Let (21) hold and let 
(31) 

u+(t) = {ut(t),ut(t), . . .,ui(t)f e and v+(t) = {vf (t) , vt (t) , . . .,v+(t)f e 



be two {possibly identical) solutions of the system (28). Then {u+{t),y'^{t))m is the 
same for each t i. e. 

m— 1 

(32) Yl ^tit)'4^ + <(t)'"m{t) = const . 

i=l 

Proof. Substituting the functions v'^ into the first m — 1 conditions and the functions 
xi+ in the last condition of (28), and then multiplying the first m — 1 conditions by u'^ 
and the last condition by f +, we get 



Subtracting the first m — 1 conditions from the last condition in the latter system, we 
get 

/ m— 1 

(33) 



f^{t)[Y,ut{t)vt{t) + uUt)v^{t)\eV' 



Since the second multiplier in (33) belongs to Vn (see (31)), (21) and (33) imply that 



m— 1 



(34) 



i=l 



We can interchange the roles of u and v in the above discussion to get in a similar 
manner that 



m— 1 



(35) 



Y.'^tm{t)-rv+{t)ul{t)eVt 



i=l 



It follows from relations (34) and (35) that the function in (32) belongs to fl Vj^, 
which implies (32). □ 

The proof of Theorem 2 (a) proceeds as follows. We search for a nontrivial polyno- 
mial solution 



(36) ^{t)^{xt{t),xt{t),...,xtit)) er 
of the system (28), where 

TV 

(37) a;+(t) = Ea,„r, i = l,2,...,m, 

n=0 



and explicitly determine the coefficients aj„. We will find such m linearly independent 
solutions of (28) (see (51) below). 

Equating all the non-positive Fourier coefficients of the functions on the left-hand 
side of (28) to zero, except the 0th coefficient of the jth function which we equate to 
1, we get the following system of algebraic equations in the block matrix form which 
we denote by S^: 



(38) 



< Tj ■ — D ■ Xj — 1, 



Tm-l ■ — D ■ X, 



m—1 



0, 



Fl • + r2 • ^2 + . . . + Fm-l ■ Xm-1 + • X„ = . 



Here the following matrix notation is used: 
(39) 





/ do 


di 


d2 ■ 


■ d^-i 


dN \ 




/ 7io 


7ii 


7i2 • • 


• 7j,Ar-i 


7iAf\ 







do 


di ■ 


• dN_2 


d^-i 




7ii 


7i2 


7i3 • • 















do ■ 


■ d^-s 


dN~2 




7i2 


7i3 


7i4 • • 















• 





do ) 




[jiN 





•• 





0/ 
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where 



N 



N 



t{z) = Y,dnz'' and Qit) = J] 7mr", ^ = 1, 2, . . . , m - 1, 



n=0 



n=-N 



(40) = (0, 0, . . . , 0)^ e C^+\ and 1 = (1, 0, 0, ... , 0)^ G C^+^ 
The column vectors 

(41) = (aicOji, ■ ■ ■ ,aiiv)^, i = l,2,...,m, 

(see (37)) are unknowns. 

Remark 1. We recall that if (Xi,X2, . . . ,X„j) defined by (41) a solution of the 
system (38), then the vector function (36) defined by (37) a solution of the system 
(28). 

We need to show that the system E>j (see (38)) has a solution for each j = 1, 2, . . . , m. 
Since f^{0) > (see (21)), -p can be represented as a power series in the neighbor- 
hood of 



n=0 



where bo = {f^{0)) ^>0, and the inverse of the matrix D is 



(42) 



(bo 


bi 


62 ■ 


■ 


bN \ 





bo 


bi ■ 


■ bN-2 


bN-i 








bo ■ 


■ bpfs 


bN-2 


Vo 





• 





bo J 



Determining X^, i — 1,2, . . . ,m — 1, from the first m — 1 equations of (38), 
(43) Xi = D^riX;^-S,jD^l, i = l,2,...,m-l, 

and then substituting them into the last equation of (38), we get 



Ti Ti Xm + r2D-^r2Xm + ... + T^-l r^_i Xm + DXm^ Tj L*"! 1 

(it is assumed that = D, i.e. the right-hand side is equal to 1 when j — m) or, 
equivalently. 



(44) (91 Gi + 92 62 + . . . + e^_i e„_i + X^ = D-^ T,- L>-i 1, 
where 

(45) Qi^D-^V,, i = l,2,...,m-l. 

For each j = 1,2, ... ,m, (44) is a hnear algebraic system of -|- 1 equations with 
(A^" -|- 1) unknowns. 
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The matrices Qi, i — 1,2, ... ,111 — 1, are symmetric since their entries are (see (45), 
(42), and (39)) 



(46) 



ei[k,i]^ei[i,k] 



foi k + l> N, 

N-{k+l) 

^ bn%,k+i+n for k + l< N. 

n=0 



Therefore ©jOj = ©i©*, i — 1,2, ... , m—1, are non-negative definite and the coefficient 
matrix of the system (44) 



(47) 



A = OiOi* + 0202* + . . . + 0„_iO;,_i + 



(which is the same for each j = 1,2, ... ,m) is positive definite (with all eigenvalues 
larger than or equal to 1). Consequently, A is nonsingular, det A > 1, and the system 
(44) has a unique solution for each j. Furthermore, A has a displacement structure of 
rank m (see Appendix B) which reduces the computational burden for solution of the 
system (44) from 0{N^) to 0{mN'^) (see [14; App. F]). 

Finding the matrix vector from (44) and then determining Xi,X2, . . . ,X^_i 
from (43), we get the unique solution of §j. To indicate its dependence on j, we 
denote the solution of by (X^, X^, . . . , Xl^_-^, X^^, 

(48) Xi :^{a^,a^,...,4j^, i = l,2,...,m, 

so that if we construct a matrix function V{t), 





vtiit) 




• ^'ltm-l(^) 




(49) V(t) = 


^m-l,l(^) 




' '^m-l,m-l(^) 














by letting (see (48)) 




N 
n=0 






(50) 




1 < i, j < m 





(note that (49) has the structure required in Theorem 2 (a); see (15), (23)), then its 
modified columns v'^{t),v'^{t), . . .,V^{t), and V^{t), 

(51) V^{t) = {vtj{t),vtj{t),...,vi_,/t),vi^{t)), j = l,2,...,m, 

are solutions of the system (28) (see Remark 1). Hence, because of the last equation 
in (28), 

(52) 

and, by virtue of Lemma 3, 
(53) {V'{t),V^{t)) 



F{t)V{t) e P+ 

- Cjj, i, j — 1,2, . . . , m. 
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for each t e T. Besides, we have 

(54) det V (t) ^ const, teT. 
Indeed, the inclusion 

(55) det V{t) e V 

is obvious (see (49) and (50)). The relation (52) implies that (see (18)) 

/+(t) detl^(t) = detF(t) detV{t) G V+. 
Thus, it follows from (21) and (55) that 

(56) detV{t)eV+. 
Next we have (see (98) below) 







(l 


■■■ 


-Ci//+ \ 









1 ■■■ 


-C2//+ 


(57) 







••• 


-C3//+ 









••• 


1 -Cm-l/f^ 








••• 


/ 


and 










(58) 


det {F- 


ya) = 





Since the column vectors (51) are solutions of the system (28), we have 



(59) 0+(t) = 0(tX(t) - f^{t)v^{t) e P+, l<j<m,l<i<m. 

Direct computations give (see (57), (49), and (59)) 



{F-y{t)v{t) = {f^{t)) 



-1 



\ t;+i(t) <2(^) 
Thus, there exists a matrix function $^(t) G (hence 
(60) det*+(^)eP+) 
such that 



{F-'y{t)v{t) = {fm))-'-^^t). 



Consequently (see (58)), 
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SO that 

(/+(t))™"'det\/(t) =det$+(t). 

Thus, it foUows from (21), (55), and (60) that 

(61) dctV{t)eV^. 

The relations (56) and (61) imply detV(t) G fl V~ yielding (54). 

The matrix function V{t) is not yet unitary, but it can be easily made such by mul- 
tiplying from the right by a constant matrix. Namely, the matrix C = {cij)ij=i^2,...,m 
defined by (53), 

(62) C = {V*{t)V{t)f, t eT, 

is nonsingular. Indeed, if C were singular and 7^ w = {wi, W2, ■ ■ ■ , Wm) £ C"* were 
such that wC = 0, then 

m 

||5]^,V^(t)||^^ = wCw* = 

i=i 

for each t G T, i.e. the vector functions V^{t), ^^(t), . . . , V™(t) would be linearly de- 
pendent. But this is impossible since the linear functional L : x x . . . x C"* 
which maps (x^ (i) , (t) , . . . ,a;+(t)) into the 0th Fourier coefficients of the functions 
standing on the left-hand side of the system (28), i.e. into (co{Ci(^)^m(^)~/^(^)^i^(^)}' 

. . ., Co{Crn-i{t)xUt) - /+(t):r+_i(t)}, Co{Ut)xt{t) + Ut)4{t) + ■■■ + 
transforms m vector functions V^{t), V"^(t), . . ., V'^{t) into linearly independent stan- 
dard bases of C", namely, \-i(y^{t)) = {6ji, 5j2, ■ ■ ■ , Sj^), j = 1,2, ... ,m, because of 
(38). Consequently, V{1) is also nonsingular since = V*{1)V{1) (see (62)). Let 

(63) U{t) = V{t){V{l))-\ 
Then U{t) is unitary since (see (63), (62)) 

U*{t)U{t) = {{v{i))-'Yv*{t)v{t){v{i))-' = {{v{i))-Yv*{i)v{i){v{i))-' =im. 

Since the matrix (^(1)) ^ is constant, U{t) G i4i(T) has the same structure (15), (23) 
as V{t), and 

(64) F{t)U{t) G P+ 

holds as well (see (52) and (63)). Moreover, detf/(t) = const, t G T (sec (54) and 
(63)), which implies that dett/(t) = 1 as we have U{1) = 1^ (see (63)). Consequently, 

(65) U{t) G 5W^(T). 
Let now 

(66) Upit) = U{t) ■ {FU{0))~^y/FU{0){FU{0))*. 

The multiplier of U{t) in (66) is a (constant) unitary matrix, so that Up{t) G Wm(T), 
it has the structure (15), (23) (since U{t) has this structure), the inclusion (24) holds 
(see (64), (66)), and (25) is valid too (see (66)). The relation detUpit) = 1 holds 
since detUpit) — c where |c| = 1 (see (65), (66)), while c > since we know that 
< det{FUF){0) = c-detF(O) = c/+(0) (see (25), Corollary 1, and (18)) and 
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/^(O) > (see (21)). Consequently, the matrix function Upit) G SUmO^) defined by 
(66) satisfies the requirements of Theorem 2 (a) and it has been constructed exphcitly. 
The proof of the part (a) is finished. 

Remark 2. Note that, as in the case ofU{t), the modified column vectors of Upit) 
are solutions of the system (28) since this property of a matrix function is preserved 
when we multiply it by a constant matrix from the right. 



6. Description of the Method 

A brief outhne of the method is the following: S{t) is approximated by MN(t)M'^{t), 
where M^it) is a lower triangular matrix function with analytic entries on the diagonal 
and whose entries below the diagonal have only finite number of nonzero Fourier 
coefficients with negative indices, the last product is represented as M^{t){M^)*{t), 
where an analytic matrix function [t) is constructed explicitly, and its convergence 
to S^{t) is proved. 

Given a matrix spectral density (1), first the lower-upper triangular factorization of 
S{t) is performed, 

(67) S{t)=M{t){M{t))\ 
where 

(68) 



M{t) 



( ftit) 

6i(t) 



Cr2{t) 







Cr,r-l{t) 










The functions fmi't)^ ni = 1,2, . . . ,r, on the diagonal are taken the canonical spectral 
factors of the positive functions det Sm{t)/ detS'm_i(t), where 5'o(i) = 1 and Sm{t) — 
{S{t))^^^, the upper-left m x m submatrix of S{t). Namely, 



(69) 

where (see (4)) 



{det S^yjz) 
(det5„_i)+(^)' 



{dot Sm^iz) 



2tt 



+ Z 



- z 



\0gdQtSra[e''^)de 



We have logdetS'm(t) G -^^i(T), m = 0,1,..., r, by virtue of (2) (see, e.g., [4; Sect. 
5]), so that the functions {deX,Sm)'^{z), and consequently f^i^), rn = 1,2, ... ,r, are 
well defined in (69). The entries ^ij, 2<i<r,l<j<i, can be found in a standard 
algebraic way from the relation (67). 



Note that (67) implies \f^\^ = Su e Li and X)!=i l^iP + I/, 



e Li, i 



2, 3, . . . , r. Thus M{t) G L2(T) (and hence M(t)U{t) G L2(T) for any U(t) G Ur(T)). 
Furthermore, G O2, m— 1,2, ... ,r, which implies that 

(70) det M{t) = ft{t)f^{t) . . . f+{t) G O^/, . 
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We search for U{t) G SUriT) such that M{t)U{t) is a spectral factor of S{t) and 
continue the description of our method in terms of Propositions 1 and 2, which can 
be easily proved using Theorems 1 and 2, respectively. 

Proposition 1. A spectral factor of S(t) can be represented as 

(71) S+{t) = M{t)U2{t)Us{t) . ..IJrit), 

where U^(t) e SUr{T) has the block matrix form 





(72) 



Urjt) 




, m = 2,3,...r-l, U,(t) = Up^it), 



Fm in (72) is the matrix function of the form (13) whose last row coincides with the 
last row of (M^_i(t))^^^, 

(73) Mi(t) := M(t), M„(t) := M {t)\J ^{t)^ ^{t) . . . U„(t) = M„,_i(t)U„,(t), 

and Up^^it) G 5Wm(T) is t/ie corresponding matrix function determined according to 
Theorem 1, m = 2, 3, . . . r. 

Proof. Obviously, the product of two matrix functions from SlAriT) is in the same 
class. Thus, by virtue of Lemma 2 (see (67), (70)), it suffices to show that 

(74) S+{t) = Mr{t) = M{t)V2it)V^{t) . . . Vrit) G L+(T). 
It follows from the structures of the matrices in (72) and (73) that 

(75) (^-W)™x™ = iMrn-l{t))^^.^UpJt), 

while the last r — m columns of M{t) remains unaltered in M^{t). 
We show by induction that 



(76) 



G Lt(T 



m 



1,2, 



Indeed, clearly (76) is correct for m 
replaced by m — 1 in it, i.e. 



1. Assume now that (76) holds when m is 



(77) (^-i(^))(.-i)x(„.-i)ei^2(T) 
Then (M^_i(t))^^^ G L^iT) has the form 



(78) (M^_i(t)), 



/^^2(^) 



At) 



V Cl(^) C2{t) 



A*m-l,m-l(^) 
Cm-lit) 



\ 






mt)j 



where iJ,±{t) G L+(T), i, j = 1, 2, . . . , m- 1, by (77), Q{t) G L2(T), j = 1, 2, . . . , m- 1, 
/+ G is defined by (69) (see (68)), and {Ci(t) , C2(t) , ■ ■ ■ Cm-i(t) , f+(t)) is the last 
row of Fm{t) according to its definition in Proposition 1 (note that Theorem 1 can 
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be applied to the matrix function Fjn{t)). The direct computation shows that (78) is 
equal to (see (13)) 

V ••• 1/ 

and we get that (see (75)) 
(79) 

{MUt))^^^ = (M.-i(t))„_f/^Jt) = (^(^--(^)V-i)x(^-i) 0^ F^m.Jt) 

belongs to C^) (^^^ (I'''))- Thus (76) is valid and taking m = r in (76) we 

get (74). 

Proposition 1 is proved. □ 

Remark 3. One can see from the above proof that (-^m(^))^x^ is a spectral factor 
of S^it), 

S^{t) = {Mmit))^^^, m = l,2,...,r. 

Thus the representation (71) realizes the step-by-step factorization of the upper-left 
submatrices of S{t). 

Remark 4. There is an alternative way of representation (71) which avoids pre- 
liminary computation of entries ^ij{t) in (68). Namely, we can compute only /^(O; 
m = 1.2, ... ,r, in (68) according to (69) {note that (Mi(t))^^^ = fi{t) ) and deter- 
mine i^^mit))^^^ recurrently from (^m-i(^)) (•^_^-)x(m-i) the formula (79). The 
entries Ci(^)) C2(^)) " " " Cm-i(^) of Fmit) and of (78) can be determined from the equation 

(-/^m-l(^))(^_l)x(m-l)(^l' C2 ■ ■ ■ , Cm-l)* = (^Im, S2m: ■ ■ ■ : Sm-l.m)^ 

which follows from (M^_i(i))^^^(M^_i(i))^^^ = Sjn{t)- In this way, we can obtain 
each {Mm{t))^^^, m = 1, 2, . . . , r, and respectively S+{t) = {Mr{t))^^^ = Mr{t). 

Relying on Proposition 1, we recurrently approximate S^{t) as follows. Let N2,N3,. . ., 
Nj. be large positive integers, and let 

(80) S+(t) = S+[N2, iVg, . . . , Nr] := M{t)lJ2{t)Vs{t) . . . U,(t), 
where Um(^) £ 5Wr(T) has the block matrix form 

(81) U„(i)- f^^-^^^^^ ° y m = 2,3,...r-l, U.(i) = C/^(..)(t), 

Fm £ -^2(1) is the matrix function of the form (13) whose last row coincides with the 
last row of (M„_i(t))^^^, 

(82) Mi(t) := M(t), M^(i) := M(i)U2(i)U3(t) ■ ■ .U^(i) = M^_i(i)U^(f), 
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Frn it) in (81) is Pj^^^^F„i{t) (sec the definition of the projection operator P/v in 
Section 3), and Up(,N^){t) G SUmiT) is the corresponding matrix function determined 
according to Theorem 2 (a), m = 2,3, ...r. We emphasize that as S{t) is given 
and the positive integers N2, N3, . . . , Nr are fixed, each Mm{t) = Mm[N2, N3, . . . , Nm], 
m = 2, 3, . . . , r, can be exphcitly constructed according to the proof of Theorem 2 (a). 
The following proposition shows that 

(83) S+(t) = S+[N2, TVs, ■ ■ ■ , Nr] = Mr(t) 
(see (80) and (82)) approximates a spectral factor of S{t). 
Proposition 2. \\S+{t) - S+{t)\\L2 ^0 as N2, N3, . . . , , Nr ^ 00. 
Proof. We prove by induction that 

(84) \\Mmit) - Mrn{t)\\L, ^ as A^s, ATg, . . . , A^^ ^ oo, m = 2, 3, . . . , r. 
Indeed, assume that 

(85) \\Mm-lit) - Mm-lit)\\L, ^0 as N2,N3,...,Nra-l^ OO 

holds (note that Mi{t) = Mi{t)). Then, by virtue of the definitions of F^{t) and 

Fm (^) ) 

(86) ||F„(t) - F^{t)\\L, ^ as N2, N3, . . . , N^-i ^ 00. 
Obviously (see (20)), 

(87) ||F„(t) - Fj^\t)\\L, ^ as iV ^ 00. 
It follows from (86) and (87) that 

\\F^{t) - F^-\t)U, ^ as N2,Ns,...,N^^oo. 

Thus, by virtue of Theorem 2 (b), Uf,(Nm){t) =r Up^it) and hence (see (72) and (81)) 

U,„(t) ^ Um(t) as N2,N3,...,Nm 00. Consequently (see (73), (82), (85), and 
Statement 1 in Sect 3), 

\\Mm{t) - Mm{t)\\L, = \\Mm-l{t)\Jm(t) - M^_i(i)U^(i) ^ 

as A^2, -/Vs, • • • , Nm oo and (84) holds. 

If we substitute m = r into (84), we get the proposition (see (74) and (83)). □ 

Remark 5. The rate of convergence in Proposition 2 is estimated under minor re- 
strictions on S{t), which is the subject of a forthcoming paper. 

Remark 6. In actual computations of S^(t) according to (80), we cannot take M{t) 
exactly since it requires scalar spectral factorizations. As it was mentioned above, our 
method does not contain any improvement in approximate computation of M{t) . We 
can assume that it can be constructed Mi{t) — Mi[Ni\{t) in (82) such that \\M{t) — 
^i(^)||l2 ^ as Ni ^ oo, and the rest of the proof of Proposition 2 goes through 
without any change. 

If we wish to construct an approximation to the canonical spectral factor 5*^, then 
we take (see Lemma 1) S+{z) = ^+(^) (^+(0))"^-^^+(0)(^+(0))* . 
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7. Numerical Simulations 

The computer code for the factorization of polynomial matrix functions by our 
method was written in MatLab in order to test the algorithm numerically and com- 
pare it with other existing software implementations available in the MatLab toolbox 
"Polyx" . The results of our numerical simulations are presented in this section. 

Two different commands, spf (•) and spf (■, syl) are available in Polyx to perform poly- 
nomial matrix spectral factorization for a discrete time variable z. (As it is explained 
in the software manual these factorizations are based on the Newton- Raphson iteration 
and on the Sylvester's method, respectively.) We have supplied the three programs 
with the same data and compared their performances. The computer with character- 
istics Intel(R) Core(TM) Quad CPU, Q6600 2.40GHz, 2.40 GHz, RAM 2.00Gb was 
used for these simulations. 

In the first place we took a test matrix whose spectral factorization was known 
beforehand, 

2^-1 + 6 + 2z 7z-^ + 22 + llz\ _ f 2 + 1 \ f 2 + z 7 + dz'" 
nz-^ + 22 + 7z 38z-^ + 8A + 38z) ~ \7 + 5z-^ 3 + \^ 1 3 + z 

(the matrix is very simple, but its determinant, —z~^ + 2 — z^, has two double zeros 
on the boundary, which usually causes difficulties in many methods). So the correct 
answer for the (right) spectral factor (with the uniqueness restriction for the coefficient 
matrix of the highest degree of z to be upper triangular with positive entries on the 
diagonal, as it is in Polyx) is 

2.000... 7.000...A ^ A.OOO... 5.000... 



1.000... 3.000... y ' \^ 1.000. 

The resulting coefficient matrices obtained by spf (•) and spf(-, syl) were the same 

'1.998657938438840 6.995302784535943\ A .000671030780579 5.001537986446627\ 
^1.002013092341746 3.007856444481516 J \ 1.001621242570816 J 

and the time elapsed varied within 0.22-0.24 sec. Below we present the results of 
computation by the program based on our algorithm which shows the advantage of 
the proposed method. In the process of the calculations three different pairs of tuning 

parameters were used: e, the accuracy level of scalar spectral factorizations of Si{t) 
and dct 5*2 (t) in (69), and N = N2, a. positive integer in (83). Accuracy improvements 
are evident as proved theoretically in Section 6: 
e = 0.0001; N = 20; time elapsed: 0.04 sec. 

'1.999540036680776 6.998506647352555\ /1.000213381614498 5.000347044583584\ 



^1.000742648478974 3.002561886914702 ^ ' \^ 1.000664442251186 J 

e = 0.000001; N = 30; time elapsed: 0.14 sec 

'1.999999670218186 6.999998845763626\ /l. 000000439718603 5.000000357258912\ 

^1.000000494672526 3.000001951212981 J '\ 1.000000164890947 J 

e = 0.00000001; = 40; time elapsed: 0.31 sec 

/1.999999957974826 6.999999852911864\ 71.000000021012597 5.000000047560575\ 

1^1.000000063037774 3.000000246615687 J \ 1.000000051966992 ) 
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When data were selected at random and exact results were unknown, the mean of 
absolute values of polynomial coefficients of the error matrix (5*+) — S was taken in 
the capacity of an accuracy estimator (in general, the closeness of [S^) to S does 

not imply that is close to S^, see [13], [1], but this is the case for polynomial matrix 
functions). In the table below this mean is denoted by £. Calculation time values are 
shown, and the matrix sizes are given; say 4x10 indicates that a 4 x 4 test matrix was 
selected with (Laurent) polynomial entries of degree 10 (with coefficients from -10 to 
10). The results of calculations by spf(-) and spf(-,s7//) were almost identical. In the 
case of our algorithm, we varied the tuning parameters of the program {N2, N^, . . . ,Nr 
in (83)) so as to obtain a slightly higher accuracy than by spf (•) and spf (•, syl), while 
the advantage in time was noticeable. 





matr 
size 


time 

sec 


accur. 

e 


matr. 
size 


time 

sec 


accur. 

e 


matr 
size 


time 

sec 


accur. 

e 


matr. 
size 


time 

sec 


accur. 

e 


spfi-) 


4x10 


0.67 


10-"' 


6x15 


5.8 


10-® 


10x20 


218 


10-® 


15x20 


1949 


10-^ 


spf{-,syl) 




0.56 






4.7 


10-® 




214 


10-® 




1952 


10-^ 


New Alg. 




0.46 


10-12 




3.4 


10-^ 




65 


IQ-® 




216 


10"® 



We express our gratitude to PhD student Vakhtang Rodonaia for working out the 
software for testing our algorithm and collecting the numerical data. 



8. Conclusion 

A new algorithm of matrix spectral factorization is developed, which factorizcs any 
matrix spectral density that admits spectral factorization. The advantage of the algo- 
rithm is illustrated by the examples of numerical simulations. 

9. Appendices 

A. Convergence properties. In this appendix we continue the proof of Theorem 

2 started in Sect. 5 and prove the second part (b), which deals with convergence 
properties of the algorithm. This proof is similar to the one given in [12] for the 
two-dimensional case. 
Observe first that: 

(i) if {Up(N){t)}Neno> No C N; is any convergent almost everywhere subsequence of 
UF(N){t), I.e. if 

Up{N){t) U{t) a.e. as No 9 A^ ^ 00, 

then 

(88) F+{t) = F{t)U{t). 

Indeed, passing to the hmit in the relations (22), (24), det {F^^^Uf(n)){z) = f^^^iz) 
(see Corollary 1 and (18)), and (25), we get 

U(t) e SUm{T), F(t)U(t) e L+(T), {FU)(z) = f+(z) e 0° (see (14)), and FU{0) > 0, 

which implies (88) (see Lemma 2). 
Now it will be shown that 
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(ii) from each subsequence {Up{N){t)}Nenicn we can extract an a.e. convergent sub- 
sequence {Up(N){t)}NenocNi- 

This will finish the proof of the relation (26) by virtue of the uniqueness of the 
canonical spectral factor and the property (i). 

We say that a sequence of functions /„ G L2, n — 1,2, . . . belongs to /C, {/n}neN G 
if one can extract a convergent in L2 subsequence from /„. Recall that an operator 
K : — >• L2 is called compact if {K{hn)}nen G ^ for any bounded sequence {/injnGN, 
\hn\ < c, n ^ 1,2, .. . (see [15; §4.6]). 

To prove the property (ii), observe that Hankel's operators 

H^:L^^ L2, C e L2, and H} : ^ Lt, f G L2, 

defined by 

(89) H^{u+)^P-{CU+), u+eL^, 
and 

(90) H}iu-) = P+ifu-), u-eL;^, 

are compact operators as limits of finite-dimensional operators (see, e.g., [15; Th. 

4-6.1]). 

Fix arbitrary j < m, and let (m^^^'', . . . , u^^i ,Um^^)'^ be the jth column of 

Up(N). Since the modified columns of Up(N){t) are solutions of the system (28) (see 
Remark 2 in Sect. 5), we have 



(91) Cr<'''' - ftN)^'''' e 1 < i < m - 1, 

and 



(92) (rur> + cf'nj"" + . . . + e-vs + ^ p+ 

It follows from the compactness of the operator (89) and (19) that 



for each i = 1,2, ...,m — 1, and thus {P {f(j^^u:^}^^)}Neni G A^, because of the 
relation (92). It follows from the compactness of operator (90) and (19) that 



as weU. Hence (see (93), (94)) 

(95) 

Since /(]v)(^) ^ f^it) and f+{t) ^ for a.a. t e T (see (14)), it follows from (95) that 
{ii+^'^^}iveNi contains an almost everywhere convergent subsequence. 
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Now we will show that the same is true for {u^^^^^} 



NeNi- I < i < m - 1. 
G /C (see (19) and Statement 1) and hence {P~{C, 



Since 



+ (JV)^ 



/C, it follows from (91) that 

(96) {P~iftN) 
as well. The compactness of the operator (90) and (19) imply that 

The relations (96) and (97) imply that {fXf^u^^^^}Neni G ^ and, consequently, an 



almost everywhere convergent subsequence can be extracted from {u^^^^} 
The proof of the property (ii) is completed and thus (26) holds. 

The proof of the remaining conditions in the part (b) continues as follows. Since 
the inverse of a matrix function F of the form (13) is 

/ 



(98) 



1 


■■ 








1 








•• 








•• 


1 


Ci 




Cm-1 


/+ 


/+ 


/+ 










(F(''))-'(()=tF-i((). 



(19) implies that 
(99) 
Hence 

U^Mt) = (F(^))-'(i) • FW(i) • C/^(.)(i) ^ F-\t)F;it) 

(see (99) and (26)) and if we denote Upit) F~^{t)F+{t), then (27) holds and the 
equation in (17) follows directly. Since each matrix function in (22) has the structure 
(15), (23), the limiting matrix function Up{t) G 5Wm(T) has the structure (15), (16). 
The uniqueness of Up{t) follows from the uniqueness of the canonical spectral factor 
and the equation in (17) since F{t) is invertible. 

B. Displacement Structure. In this section we prove that the matrix A defined 
by (47) has a displacement structure of rank m with respect to i.e. (see [14; App. 
F.l]) 

(100) i?zA A - ZAZ* 

has rank m, where Z is the upper triangular (A^ + 1) x {N -\- 1) matrix with ones on the 
first up-diagonal and zeros elsewhere (i.e. a Jordan block with eigenvalue 0). There 
are several forms of displacement structure and we have selected a suitable one. 
Obviously, In+i has the displacement structure of rank 1, namely. 



(101) 



R7.I 



ZJ-N+l 



= 7 



N+l 



ZIn+iZ* = ££* 
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where E — (0, 0, . . . , 0, 1)^ e C^"*"^. We will show that for each Toephtz-hke matrix 

/ ^0 ?7i ?72 

^1 ?72 ?73 

^2 ?73 ?74 



(102) 



e 











yr/n 

the matrix 00* has the displacement structure of rank 1, namely 
(103) i?z(00*) = 00* - Z00*^* = AA*, 

where A = (770, ?7i, ■ ■ ■ , "^n)^ ■ Indeed, it follows from the definitions of matrices 
and A and from the structure of that 



Z*£: = and £:^Z = 0^, 
Z0 = 0Z* and Z0* = 0*Z*, 



(104) 
(105) 

and 

(106) - Z0Z = Al^ and 0* - Z*0*Z* = lA*, 

where and 1 are defined by (40). Since 1^1 = 1, it follows from (106) that 

(0 - ZeZ){Q* - Z*Q*Z*) = AA*. 
Hence, taking into account (105), 

AA* = 00* - ez*e*z* - zeze* + zezz*e*z* = 
00* - zee*z* - zee*z* + zezz*e*z* 

and (103) holds since (see (101), (105), and (104)) 

-zQe*z* + zQzz*e*z* = ze{zz* - In+i)Q*z* 

= -ZQ££'^Q*Z* = -QZ*££'^ZQ* = -00 0^0* = 0. 

Every matrix 0j, i — 1,2, ...,m — 1, defined by (45) has the structure (102) by 
virtue of (46). Thus we can write (103) for each i, 



(107) 



Rzi^m = ©i©* - ^©^©*^* = A^A* , i = 1, 2, . . . , m - 1. 



Since Rz defined by (100) is linear, Rz{/S.i + A2) = Rz^i + Rz^2-i it follows from 
(47), (107), and (101) that 



' m— 1 



m—1 



m—1 



i=l 



i=l 



where ^4 = [Ai, A2, . . . , h-m-i, ^\ is the (A^ + 1) x m matrix (of rank at most m) with 
columns Ai, A2, . . . , A^_i and £. 
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